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" Abstract, we describe an algorithm for the computation of generalized (or weighted) 

Ehrhart series based on Stanley decompositions as implemented in the offspring Nmzln- 
tegrate of Normaliz. The algorithmic approach includes elementary proofs of the basic 
results, we illustrate the computations by examples from combinatorial voting theory. 

>: 
o. 

Let M C Z" be an affine monoid endowed with a positive Z-grading deg. Then the 
psj ■ Ehrhart or Hilbert series is the generating function 

cm: 

, — ,: E M (t) = £ ? degx = £ #{x £ M : deg* = k}t k , 

O ■ xeM k=0 

and E(M,k) = #{x £ M : degx = k} is the Ehrhart or Hilbert function of M (see Hi for 
terminology and basic theory). It is a classical theorem that £m(0 is the power series 
expansion of a rational function of negative degree at to = and that E(M,k) is given 
by a quasipolynomial of degree rankM — 1 with constant leading coefficient equal to the 
(suitably normed) volume of the rational polytope 

> ■ P = cone(M) DAi 

t^. . where cone(M) C K is the cone generated by M and A\ is the hyperplane of degree 1 

points. In the following we assume that 

h ■ M = cone(M)HL 

t— i \ 

for a sublattice L of Z". Then E(M,k) counts the L-points in the multiple kP, and is 
therefore the Ehrhart function of P (with respect to L). 

Monoids of the type just introduced are important for applications, and in some of 
them, like those discussed in Section [31 one is naturally led to consider generalized (or 
H : weighted) Ehrhart series 



EM,f{t) = £ f{x)t^ x 

where / is a polynomial in n indeterminates. It is well-known that also the generalized 
Ehrhart series is the power series expansion of a rational function; see [Q]|, CO- 

In the last months we have implemented an offspring of Normaliz [6] called Nmzlnte- 
grateQthat computes generalized Ehrhart series. The input polynomials of Nmzlntegrate 
must have rational coefficients, and we assume that / is of this type although it is math- 
ematically irrelevant. This note describes the computation of generalized Ehrhart series 



2010 Mathematics Subject Classification. 52B20, 13F20, 14M25, 91B12. 

1 Nmzlntegrate will be uploaded to J6] together with Normaliz 2.9 by February 2013. 

1 



2 



W. BRUNS AND C. SOGER 



based on Stanley decompositions lfT3l . Apart from taking the existence of Stanley de- 
compositions as granted, we give complete and very elementary proofs of the basic facts. 
They follow exactly the implementation in Nmzlntegrate (or vice versa). 

The generalized Ehrhart function is given by a quasipolynomial q{k) of degree < 
deg/ + rankM— 1, and the coefficient of £ de g/+ rankM -! m q(j^ can easily be described as 
the integral of the highest homogeneous component of / over the polytope P. Therefore 
we have also included (and implemented) an approach to the computation of integrals of 
polynomials over rational polytopes in the spirit of the Ehrhart series computation. See 
B21 and (H for more sophisticated approaches. 

Acknowledgement. We gratefully acknowledge the support we received from John Ab- 
bott and Anna Bigatti in using CoCoALib 0, on which the multivariate polynomial al- 
gebra in Nmzlntegrate is based. 

1 . The computation of generalized Ehrhart series 

Via a Stanley decomposition and substitution the computation of generalized Ehrhart 
series can be reduced to the case in which M is a free monoid, and for free monoids one 
can split off the variables of / successively so that one ends at the classical case M = Z+. 
We take the opposite direction, starting from Z+. 

1.1. The monoid Z+. Let M = Z+. By linearity it is enough to consider the polynomials 
f(k) = k m , k E Z_|_, for which the generalized Ehrhart series is given by 

00 

£ k m t um , u = deg 1 . 
k=0 

and if necessary we can assume u = l, substituting t !->■ t u in the final result. 
The rising factorials 

(Jfe+l) m = (k+l)---(k + m) 
form a Z-basis of the polynomial ring Z[k]. Therefore we can write 

m 

J=o 

and use that 

00 00 00 / -t \ {fj J 

E(*+i),«'=E(^w = j:(OW= (—) =j T ^yn- 0.2) 

k=0 j=r j=0 V 1 V v 1 V 

Equations (11.11) and (11.21) solve our problem for M = Z + and f(k) = k" 1 : 

i^t uk = { ^l v A m , u (t)eZ[t}. (1.3) 

It is enough to compute A m ^(t) because A mtU (t) =A m ^(t u ). One should note that A m ^ u is 
a polynomial of degree m. Therefore the rational function in (11.31) has negative degree. 

Since the coefficient s m ^ m of (k+ \) m in the representation of k m is evidently equal to 1, 
we have 

00 1 

2^ k m t um = — — ^ + terms of smaller pole order at t = 1 (1 .4) 

k=o {l-t) m 
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Remark 1. The coefficients s m j in (11.11) and the coefficients of the polynomials A m i are 
well-known combinatorial numbers. 

(a) s m j = (—l) m ~jS(m+ 1,7+1) where S{p,q) is the Stirling number of the second 
kind that counts the number of partitions of a p-set into q blocks. This follows imme- 
diately from the classical identity k m+l = U]=t{-\) m+l ~ j S(m+\J)(k)j (for example, 
see Stanley 03 4.3,c]). 

(b) For m = we have Aoj = 1 and A m _\ = LT=l A(m, j)t-> for m > where A(m, j) is 
the Eulerian number lfT4l 4.3,d]. 

1.2. The monoid Z d + . Next we consider M = Z+. The crucial observation is that the 
problem is multiplicative for products of polynomials in disjoint variables. Suppose that 

/(■*) = g(y) h (z),y= (x h ...,x r ),z = (x r+ i,...,x d ). Then 

Em j{t) = £ /(*)**«* = ( £ g(y)t d ^) ( £ h(z)t de A d.5) 

by multiplication of power series. 

In order to exploit (11.51) we split the last variable off, 

f( x ) = Y,fi(xi,---,Xd-l)x l d , 

i 

and obtain 



(' 



1 r J yezf / 



(1.6) 



with w = degej. 

Applying this formula inductively allows us to eliminate all variables xi and to end with 
the desired representation of Eyd ^(t). 

Generalizing (11.41) . let us consider the case in which / is a monomial, /(jq, . . . ,x d ) = 

x™ 1 ■ ■ ■x™'*, and 7L\ is endowed with its standard degree, deg(x) = x\-\ \-x d . Then 

equations (11.51) and (11.41) imply that 

ffX\\ ' ' ' fflcl 

Em fit) = -, r — ; — — — —, + terms of smaller pole order at t = 1 . (1.7) 

' rl u \ ) q Ami-\ hm^+a r 

1 .3. Using the Stanley decomposition. We now turn to general McZ". Normaliz com- 
putes a triangulation E of cone(M) into simplicial subcones a. Moreover, it computes a 
disjoint decomposition 

cone(M) = (J C7\5 a 

where S a is a union of facets of o. The existence of such a decomposition is a nontrivial 
fact. Classically it is derived from the Bruges ser-Mani theorem on the existence of line 
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shellings (see Stanley lfT3lO . Instead of a line shelling, Normaliz (now) uses a method of 
Koppe and Verdoolaege: see IITOll and E Section 4]. 

Every simplicial subcone (of full dimension) is generated by linearly independent vec- 
tors vi , . . . , v d E M, d = rankM. They generate a free submonoid M a of M. For every a 
Normaliz computes the set 

E a = {xegp(M) :x = aivi + --- + a d v d: a t e [0,1)}. 

For x G E a we let e(x) be the sum of those v ; for which (i) a, = and (ii) the facet of o 
opposite to Vj lies in the excluded set S a (so that x lies in the excluded set). Then it is not 
hard to see that we have a disjoint decomposition 

M = U U x + £{x)+M a . 

aez xeE a 

It is called a Stanley decomposition since its existence is originally due to Stanley lfT3l . 
In the following we set x = x + e(x) and 

N ^ = x + M a . 

Then 

%W=I E E Na ,j(t). 

Set d = rankM, and for given o consider the linear map 

a a : Z c l ->Z n , a a (y h ...,y d )=yiv 1 + --- + ---+y d v d , 
where vi , . . . , v d is the generating set of M a as above. With 

deg a v = dega a (y), 
8aM=f(a„(y)+x), (1.8) 

we have 

This equation transforms the summation over N(j^ x into a summation over Z d + . Then we 
can apply (11.61) inductively to 

£a,/W= £ E Na ,j(t). (1.9) 

Finally, we sum the rational functions E a j(t) over the triangulation E. 

Remark 2. (a) Instead of applying (11.61) to every a, we accumulate the polynomials g CjX 
over all o that induce the same degree deg CT on 7L d (the classes formed in this way are 
called denominator classes). 

(b) The time critical steps in the algorithm are 

(1) the coordinate transformation (|1.8I) . and 

(2) the inductive application of (11.61) . 
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In order to speed up (1), we factor the polynomial /, transform the factors separately, 
and multiply the transformed factors. If / happens to decompose into linear factors, then 
multiplication of linear polynomials becomes a time critical step. In order to speed up (2) 
we have introduced the denominator classes. 

(c) Note that Y, ye z<t 8a,x(y)t ° y is invariant under permutations of variables yi that 
preserve the degrees deg CT e ; . Therefore one can go over g ax monomial by monomial and 
reorder the exponent vectors in such a way that the exponents of variables corresponding 
to the same degree become decreasing. The reordering significantly reduces the number 
of monomials in the polynomials to which (11.61) must be applied, saves memory and also 
speeds up (11.61) . 

(d) We want to point out that (11.61) is not applied recursively. Instead the right hand 
side is expanded after the elimination of xj, and Xd-\ is then eliminated from the resulting 
polynomial whose coefficients are rational functions in t. This procedure is repeated until 
all Xj have been eliminated. 



2. The quasipolynomial, its virtual leading coefficient, and 

integration 

2.1. The quasipolynomial. All rational functions in t that come up in (11.91) can be writ- 
ten over the denominator 

( \ _ ji'^deg/+rankM 

where t is the least common multiple of the numbers degx for the generators x of M that 
appear in the triangulation. This follows from (|1.6I) if one observes that 1 — t u divides 
1 — r. Moreover, all summands have negative degree as rational functions in t. Therefore 
|[T4l 4.4.1] implies the following proposition. 



Proposition 3. 



k=0 



where q is a rational quasipolynomial of period % dividing £ and of degree < deg/ + 
rankM- 1. 

The statement about the quasipolynomial means that there exist polynomials qV\ j = 
0, . . . , n — 1 , of degree < deg / + rankM — 1 such that 

q(k) = qU\k), j = k (*), 

and 

a <J)(k) - J J) +a U) k+- ■ -+a U) r deg/+rankM-l 

with coefficients qf' G Q. As we will see below, it is justified to call 

ed(M,/) = deg/ + rankM - 1 

the expected degree of q. 
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2.2. The virtual leading coefficient and Lebesgue integration. Let m = deg/ and 
write / = f m +g where f m is the degree m homogeneous component of m. Then degg < m, 

and it follows from Proposition [3] that g does not contribute to the coefficient %l M jy 
Moreover, this coefficient is independent of j and given by an integral, as we will see in 
Propo sition ID below. 

For the representation as an integral we must norm the measure in such a way that it is 
compatible with the lattice structure. We will integrate over the polytope 

P = cone(M)nAi, A 1 = {x G R" : degx = 1}. 

Let Lq — Lfl MM DAq where Aq = {x G W 1 : degx = 0} is the linear subspace of degree 
elements. Then Lq is a (saturated) sublattice of L of rank d — 1 (d = rankM), and 
we choose a basis u\ , . . . , u c i-i of Lq. Note that H = MM C\A\ has dimension d — 1 and 
contains a point z G L since we have required that deg takes the value 1 on gp(M), and 
we can consider the basic Lo-simplex 8 = conv(z,z + «i, • • • ,Z + w^-i) in H. Now we 
norm the Lebesgue measure A on H by giving volume l/(d — 1) ! to the basic Lo-simplex. 
(The measure is independent of the choice of 8 since two basic Lo-simplices differ by an 
affine-integral automorphism of H.) We call A the L-Lebesgue measure on H. 

Proposition 4. For all j = 0, . . . , % — 1 one has 

9%ij) = J p UdX. (2.1) 

Proof. We may assume that / is homogeneous of degree m. Let 

1 

L c = -L. 

c 



Then 



/ fmdl = lim V -inrf( x 

JP xeP<lL c c 



by elementary integration theory. 
Note that 



by homogeneity and that x G P D L c if and only G Lfl cP. Thus 

J p f m dX = \im £ ^n/(y). 



On the other hand, we obtain q^l M ^ as the limit over the subsequence (bn + j)bel+ : 



by Proposition^ This concludes the proof. □ 

In view of Proposition |4] it is justified to call q s d(M.f) = tfJ^M f) me v ^ rtua ^ leading 
coefficient, and the proposition justifies the term "expected degree" for deg / + rankM — 1 
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the. In analogy with the definition of multiplicity in commutative algebra (for example, 
see fl5j), we call 

vmult(M,/) = ed(M,/)!? ed(Af)/) 

the virtual multiplicity of (M, /) . It is an integer if P is a lattice polytope and f m has 
integral coefficients, as we will see below. 

2.3. Computing the integral. It is natural to compute the integral by summation over 
the triangulation: the triangulation of cone(M) into simplicial subcones o induces a tri- 
angulation of the polytope P into simplices 8 = o fl P. As usual let vi , . . . , vj G M be the 
generators of a. Then 8 is spanned by the degree 1 vectors v,/ deg(v,-), i = 1,. . . ,n. Let 
e\ , . . . , be the unit vectors in W l . Then the substitution e,- \-± v,-/ deg(v,) induces a linear 
map M. d — > MM that in its turn restricts to an affine map a from the standard degree 1 
hyperplane in M. d spanned by e\ , . . . , to the hyperplane H = A\ n MM, and the image of 
the unit simplex A is just 8. 



Proposition 5. One has 



fik = |de»,(v, ,v,)| l(foa)d (2 . 2) 



deg(vi)---deg(v d ) 7a 

where /i «5 Z, d -Lebesgue measure on the hyperplane H of standard degree 1 m M rf anJ 
detL(vi, . . .,Vd) is the determinant of the coefficient matrix ofv\, . . . with respect to a 
basis ofL fl MM. 

Proof. This is just the substitution rule if one observes that the absolute value of the 
functional determinant of a\H is given by the factor in front of the integral. For an affine 
map the functional determinant is constant. So we can assume / = 1 and it remains to 
relate the volumes of 8 and A. But A has volume 1 / ( d — 1 ) ! with respect to /i and 8 has 
volume 

1 |det L (vi,...,v rf )| 
(d-l)ldeg(v 1 )---deg(v d ) 
with respect to A; see B2l Section 4]. □ 

After the substitution it remains to evaluate the integral over A, and this can be done 
monomial by monomial: 



Proposition 6. 



A (mH \-m d + d-l)\ 



Proof Let g = y mi ■ ■ ■y 1 / and M = Z+. Then 

m\ ! • • • m d \ 

Em e(t) = r « + terms of smaller pole order at t = 1, 

'« w (1 — t) ( m] H — 

as stated in (11.71) . 

The quasipolynomial is a true polynomial in this case, and the (virtual) multiplicity is 
given by the value of the numerator polynomial at t = 1, namely mi ! • • • m d \ (for example, 
see (51 4.1.9]). Now Proposition|4]gives the integral. □ 
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3. Computational examples 

We illustrate the use of Nmzlntegrate by three related examples coming from combi- 
natorial voting theory that are discussed in lfT2l . We refer the reader to ffTTTl . lfT2ll or lfl"5ll 
for a more extensive treatment. 

Consider an election in which each of the k voters fixes a linear preference order of n 
candidates. In other words, voter i chooses a linear order of the candidates 1,. . . ,n. Each 
such order represents a permutation of 1 , . . . , n. Set N = n\. The result of the election is an 
AMuple (x\,.. . ,jiqv) in which x p is the number of voters that have chosen the preference 

order labeled p. Then x\-\ \-x^ = k, and (x\, . . . ,Xfj) can be considered as a lattice 

point in the positive orthant of or, more precisely, as a lattice point in the simplex 

t/W = nA k = kiml n Ai) = 

where is the hyperplane defined by x\ -\ \-x^ = k, and = U^ 1 ' is the unit 

simplex of dimension N — 1 naturally embedded in Af-space. We assume that all lattice 

(n) 

points in the simplex U} have equal probability of being the outcome of the election. 
The following three problems have been considered in lfl2l for 4 candidates A, 5, C,D: 

(1) the Condorcet paradox, 

(2) the Condorcet efficiency of plurality voting, 

(3) plurality voting versus cutoff. 

For n = 4 one has N = 24, and the dimension of the polytope is already quite large. 

Let us say that candidate A beats candidate B if the number of voters that prefer can- 
didate A to candidate B is larger than the number of voters with the opposite preference. 
Candidate A is the Condorcet winner if A beats all other candidates. As the Marquis de 
Condorcet noticed, the relation "beats" is nontransitive for some outcomes of the election, 
and there may be no Condorcet winner. This phenomenon is called the Condorcet para- 
dox. Problem (1) asks for its asymptotic probability as the number k of voters goes to oo, 
or even for the precise number of election results without a Condorcet winner, depending 
on the number k of voters. 

It is not hard to see that the outcomes that have A is the Condorcet winner can be 
described by three homogeneous linear inequalities A/(x) > whose coefficients are given 
in Table Q] (relative to the lexicographic order of the permutations of A,B,C,D). They 

Xr. l l l l l 1-1-1-1-1-1-1 l l-i-i l-i l l-i-i l-i 
A 2 : 1 1 1 1 1 1 1 1-1-1 1-1-1-1-1-1-1-1 1 1 1-1-1-1 
X 3 : 1 1 1 1 1 1 1 1 1-1-1-1 1 1 1-1-1-1-1-1-1-1-1-1 

Table 1 . Inequalities expressing that A beats the other 3 candidates 

cut out a rational polytope from £A B ), and the probability of Condorcet's paradox can 
be computed from the volume of the polytope. Finding the precise number of election 
results without (or with) a Condorcet winner requires the computation of the Ehrhart 
function of the semi-open polytope . Neither Normaliz nor Nmzlntegrate can yet compute 



THE COMPUTATION OF GENERALIZED EHRHART SERIES IN NORMALIZ 9 

Ehrhart series for semi-open polytopes directly, but it is always possible to fall back on 
inclu sion/exclu sion . 

We refer the reader to Q for a description of problems (2) and (3) and for the systems 
of linear inequalities to be solved in each case. Normaliz 2.8 can indeed compute the 
volumes and the Ehrhart series in dimension 24 that arise from tasks (1), (2) and (3) 
despite the fact that the triangulations to be evaluated for (2) and (3) are formidable (see 
Table[3]or Q). 

As Schiirmann lfl2l observed, the computations can be considerably simplified by ex- 
ploiting the symmetries in the inequalities: some variables share the same coefficients 
in each inequality, for example the first 6 variables in Table [TJ Therefore they can be 
replaced by their sum, and the replacement constitutes a projection of the original poly- 
topes, monoids or cones onto objects of smaller dimension. For the Condorcet paradox 
the system of inequalities reduces to Tabled However, instead of simply counting lattice 

1-1 1 1 1-1-1-1 
1 1-1 1-1 1-1-1 
1 1 1-1-1-1 1-1 

Table 2. Inequalities exploiting the symmetries in Tabled] 



points, one must now count them with their numbers of preimages. These are given by 
polynomials, namely products of binomial coefficients. In our example the polynomial is 



yi+5 
5 



(yi + 1) (y 3 + 1) (j4 + 1) (ys + 1) (ye + 1) (y? + 1; 



ys + 5 
5 



where y\ =x\-\ h X(, etc. In other words, the Ehrhart function (or the volume) of a high 

dimensional polytope is replaced by a generalized Ehrhart function of a polytope of much 
lower dimension (or the virtual leading coefficient of the quasipolynomial). 

A priori it may not be clear that the replacement of combinatorial complexity in high 
dimension by multivariate polynomial arithmetic in low dimension pays dividends, but 
this is indeed the case. Tables [3] and |4] compare both approaches. The computations were 
run on a SUN xFire 4450 with 20 parallel threads. If the computations in Table [3] are 
restricted to volumes, they become faster by a factor of approximately 3. 



computation 


triangulation size 


real time 


Condorcet paradox 


1,473,107 


00:00:30 h 


Condorcet efficiency 


347,225,775,338 


218:13:55 h 


plurality vs. cutoff 


257,744,341,008 


175:11:26 h 



Table 3. Computation times (real) for Ehrhart series in dimension 24 



A welcome side effect of the computations of the generalized Ehrhart functions is that 
they have confirmed the results reported on in ||7]|. 
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computation 


rank 


deg/ 


triangula- 


Normaliz 


gen Ehrhart 


lead coeff 








tion size 


time 


series time 


time 


Condorcet paradox 


8 


16 


17 


0.05 sec 


5.2 sec 


0.08 sec 


Condorcet efficiency 


13 


11 


17,953 


0.41 sec 


5:49:29 h 


1:54:35 h 


plurality vs. cutoff 


6 


18 


3 


0.06 sec 


18.4 sec 


0.54 sec 



Table 4. Computation times (real) for symmetrized data 
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